#!/usr/bin/bash
echo "usage: vep.sh inputfile outputdir ref"

set -o nounset -o errexit -o pipefail

if [ $# -ne 3 ]; then
   echo "!!!usage: vep.sh inputfile outputdir ref"
   exit
fi


inputfile=$1
outputdir=$2
ref=$3
vep=/picb/lilab/tools/bin/vep #公共VEP路径
name=`basename -s .vcf.gz $inputfile` #提取input的basename,如果输入文件为其他后缀，则将-s后面的参数改为后缀名

$vep \
 --species homo_sapiens \
 --assembly GRCh38 \
 --input_file $inputfile \
 --format vcf \
 --output_file ${outputdir}/${name}.vep.vcf.gz \
 --cache --dir_cache /picb/lilab/annotation/VEP/cache_101 \
 --fasta $ref \
 --distance 5000 \
 --cache_version 101 \
 --everything --allele_number --vcf --offline --force_overwrite \
 --plugin LoF,loftee_path:/picb/lilab/annotation/VEP/LOFTEE/loftee-grch38,human_ancestor_fa:/picb/lilab/annotation/VEP/LOFTEE/human_ancestor_grch38.fa.gz,conservation_file:/picb/lilab/annotation/VEP/LOFTEE/loftee.sql,gerp_bigwig:/picb/lilab/annotation/VEP/LOFTEE/gerp_conservation_scores.homo_sapiens.GRCh38.bw \
 --dir_plugins /picb/lilab/annotation/VEP/LOFTEE/loftee-grch38 \
 --plugin CADD,/picb/lilab/annotation/CADD/data/prescored/GRCh38_v1.6/no_anno/whole_genome_SNVs.tsv.gz,/picb/lilab/annotation/CADD/data/prescored/GRCh38_v1.6/no_anno/gnomad.genomes.r3.0.indel.tsv.gz \
 --custom /picb/lilab/annotation/gnomad/gnomAD/r3.x/gnomad.genomes.r3.0.sites.vcf.gz,gnomADv3.0,vcf,exact,0,AF,AF_eas \
 --custom /picb/lilab/annotation/clinvar/GRCh38/clinvar_20210718.vcf.gz,ClinVar210718,vcf,exact,0,CLNSIG,CLNREVSTAT,CLNDN \
 --plugin SpliceAI,snv=/picb/lilab/annotation/SpliceAI/spliceai_scores.raw.snv.hg38.vcf.gz,indel=/picb/lilab/annotation/SpliceAI/spliceai_scores.raw.indel.hg38.vcf.gz \
 --compress_output bgzip

##将vep的注释结果整合成txt格式
python2=/picb/lilab2/anaconda3_1/envs/python27/bin/python
parse_vep=read_vep_vcf.py 
${python2} ${parse_vep} --vcf ${outputdir}/${name}.vep.vcf.gz > ${outputdir}/${name}_vep.txt